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Ergodicity is proved for granular contact forces. To obtain this proof from first principles, this 
paper generalizes Boltzmann’s stosszahlansatz (molecular chaos) so that it maintains the necessary 
correlations and symmetries of granular packing ensembles. Then it formally counts granular contact 
force states and thereby defines the proper analog of Boltzmann’s H functional. This functional is 
used to prove that (essentially) all static granular packings must exist at maximum entropy with 
respect to their contact forces. Therefore, the propagation of granular contact forces through a 
packing is a truly ergodic process in the Boltzmannian sense, or better, it is self-ergodic . Self- 
ergodicity refers to the non-dynamic, internal relationships that exist between the layer-by-layer 
and column-by-column subspaces contained within the phase space locus of any particular granular 
packing microstate. The generalized H Theorem also produces a recursion equation that may 
be solved numerically to obtain the density of single particle states and hence the distribution 
of granular contact forces corresponding to the condition of self-ergodicity. The predictions of 
the theory are overwhelmingly validated by comparison to empirical data from discrete element 
modeling. 
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I. INTRODUCTION 

It is interesting that the distribution of granular con- 
tact forces, Pf(f ), appears to have an asymptotically ex- 
ponential tail [1-3]. In statistical mechanics, an expo- 
nential tail appears in the distribution of positive-valued 
variables (like particle energies) when (1) the system- 
wide sum of those variables is conserved, and (2) the 
ensemble has a flat measure. In thermal systems like a di- 
lute, classical gas, the first of these conditions is automat- 
ically satisfied if the particle dynamics conserve energy 
one collision at a time. The second condition requiring a 
flat measure is more difficult to prove. There is an ergodic 
theorem based on the Poincare cycle, which predicts the 
flat measure in that all accessible states will be visited 
with equal probability. Unfortunately this Poincare cycle 
has a period greater than the age of the universe and says 
nothing about the measure over timescales that we ob- 
serve. Boltzmann’s H Theorem, on the other hand, pre- 
dicts the flatness in the relevant time scale, but it does so 
at the cost of introducing a new a priori assumption, al- 
beit a very reasonable one: that colliding particles are not 
correlated before they collide. Therefore, the two-point 
correlation function describing two particles colliding is 
written as 

F(pi,P2,t) = f(p u t) ■ f{p2,t) (1) 

where / is the single particle density of states. This 
is Boltzmann’s stosszahlansatz , or his “assumption of 
molecular chaos.” The subsequent proof found in his 
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//-Theorem (not repeated here) tells us that an ensem- 
ble of systems describable by Eq. 1 after following even 
a very tiny segment of the Poincare cycle will have a dis- 
tribution of momenta (and energies) that is the same as 
the average over the entire Poincare cycle. Therefore, if 
Boltzmann’s stosszahlansatz is valid, the second condi- 
tion listed above can be relaxed: we do not need a flat 
measure over the entire accessible region, we need only a 
little bit of dynamics on any part of the accessible region 
to obtain the necessary flatness in the practical sense. 

So if this is what causes the exponential tail in a ther- 
modynamic system, then what causes the tail in a static 
granular packing? As it turns out, it is caused by the 
same two conditions listed above. This paper shall ex- 
plain the unique character of the “conservation laws” in 
a granular packing and how they are analogous to energy 
conservation in a thermal system, but with an important 
difference. This paper shall also explain what part of the 
physics causes the contact forces in a granular ensemble 
to spread out with a flat measure. Following Boltzmann’s 
lead, we do not really need to assume the measure is flat 
over the entire accessible part of phase space. All we 
need to do is write a transport equation that explains 
how the stress state of individual grains connects to the 
stress state of their neighbors, and develop an analogous 
version of Boltzmann’s stosszahlansatz and H Theorem. 

The goal in this paper shall be to explain the ideas 
physically, not just mathematically, and with as much 
clarity as possible. To do this, Sec. II will look at several 
models with increasing complexity to identify the essen- 
tial parts of the physics that determine the statistical 
mechanics of granular contact forces. After we have seen 
them, in Sec. Ill it will be a simple matter to write the 
equations and produce the generalized stosszahlansatz. 
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It will be asserted as a very reasonable assumption from 
first principles, just like Boltzmann’s version, but it will 
also be supported by empirical observations. Sec. IV will 
then write the transport equation for single grain states. 
Assuming the generalized stosszahlansatz is valid, it will 
formally prove that essentially every accessible packing 
state must exist at maximum contact force entropy, and 
thus each packing state has the same statistical proper- 
ties as an ergodic average over all states in the ensem- 
ble. This will be analogous to Boltzmann’s H Theorem, 
although it will be adapted to the unique features of a 
granular packing. Thus, we will have achieved an ergodic 
proof in the Boltzmannian sense. However, it is better 
not to call it ergodicity , but rather self-ergodicity because 
a static granular packing, being static , does not explore 
phase space. The ergodicity consists of a special set of re- 
lationships that exist among the values of its phase space 
coordinates. The proof of this self-ergodicity stands or 
falls with the validity of the stosszahlansatz. There- 
fore, Sec. V will test the predictions empirically through 
comparison with discrete element modeling (DEM) data. 
The evidence overwhelmingly validates the theory. This 
paper leads to the immediate conclusion that (in the ther- 
modynamic limit) every set of first shell states is equally 
probable in a granular packing, subject to the conserva- 
tion laws. This opens the door to wider applications of 
statistical mechanics in granular media. 


II. THE UNIQUENESS OF GRANULAR 
CONTACT FORCE PHYSICS 

A. The q Model: Ergodicity in One Dimension 

A granular packing is static and cannot explore its 
phase space. So how then does the distribution of con- 
tact forces statistically relax to have an asymptotically 
exponential tail? The q model [4] begins answering this 
question. A particular example of a q model lattice is 
illustrated in Fig. 1, consisting of a regular lattice (di- 
amond in this case) with grains at the nodes and grain 
contacts along the legs. The vertical load borne by the 
grains in the i th layer are randomly redistributed to the 
grains in the ( i 4- l) st layer by means of the random q^ 

variables , 0 < q^ < 1, which tell what fraction of the 
vertical load goes into the contact to the grain below it 
to the right, the remainder going to the grain below it 
to the left. After some 20 or so layers the model begins 
to relax so that the distribution of vertical forces has an 
asymptotically exponential tail. In this model, the ex- 
ploration of phase space is not through the time dimen- 
sion, but through a spatial dimension. All the vertical 
forces in one layer of the lattice form the coordinates in 
a single-layer phase space, and translating through the 
lattice layer-by-layer does explore that space. 

It may be noted that the symmetry of the q model may 
be improved by considering the case without gravity so 



FIG. 1: A two-dimensional diamond lattice version of the q 
model. Forces propagate layer by layer and are redistributed 
in the next layer by the q values, which play the role of the 
collision matrix. This forms an analogy to the molecular ki- 
netics of a gas in which the vertical dimension represents time 
and the complete set of molecules is a single horizontal layer. 


there is no additional weight wo applied to each grain 
(the injection term). Then, after a q model is developed 
over many layers, the upper, unrelaxed levels may be dis- 
carded, keeping only a block of deeper layers that are sta- 
tistically relaxed. Then it is impossible to tell which way 
the information was propagated across the lattice: up or 
down become indistinguishable. Indeed, one can apply q 
values and propagate forces in one of these two directions 
and then when the model is relaxed extract the apparent 
q values as if the forces had been propagated in the op- 
posite direction. Call these the p values. The author has 
incidentally noted that for continuous distributions of the 
q variable rj(q) the distribution of the p values automat- 
ically relaxes to i ?(p), identical to the distribution of the 
q values. Therefore, the weightless q model does become 
truly symmetric under inversion when it is relaxed in all 
its layers. This is analogous to a box of grains in weight- 
lessness, squeezed on all sides by the container walls. An 
ensemble of such packings must be symmetric under in- 
version (plus parity reversal in the general, sheared case), 
and hence are pictured well by the relaxed, weightless q 
model. 

Coppersmith, et a/, showed that for the case of the 
uniform q model where 77 (q) = 1 for 0 < q < 1, the 
exploration of the vertical forces’ phase space is truly er- 
godic in that it visits every state with equal likelihood. 
On the other hand, the non-uniform q models are defi- 
nitely not ergodic. Since real granular packings produce 
a distribution of Cartesian contact forces unlike that of 
the uniform q model, indicating non-uniform r)(q ), how 
could we claim ergodicity for a real granular packing? 
The answer is that the scalar q model allows too much 
freedom because it does not treat the full vectorial forces, 
and we shall see below that in the fully vectorial lattice 
of a granular packing an ergodic exploration of vectorial 
phase space is not identical to each Cartesian component 
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ergodically exploring its own phase space independently. 
But thus far in the argument, ergodicity has not been 
proven for real granular packings, neither for their Carte- 
sian components of force nor for their forces vectors. We 
shall say that the q model is a picture of randomization 
in a granular packing rather than ergodicity until we can 
prove otherwise. 


B. Analogy to Thermodynamics 

Coppersmith, et a/, introduced an analogy between the 
q model and the molecular kinetics of a classical dilute 
gas to further illuminate the nature of this layer-by-layer 
randomization. The analogy compares the static system 
of contact forces to the dynamics of a dilute classical gas: 
the vertical dimension of the q model’s lattice is analo- 
gous to the time dimension, and a single horizontal layer 
in the lattice (a subspace of the entire lattice’s density 
of states) is considered to be the entirety of the ther- 
modynamic system at any given moment of time. The 
“atoms” are the contacts between grains rather than the 
grains themselves. The “collisions” between atoms are 
the grains rather than the contacts. That is, a grain 
is a collision of the forces that exist on its contacts. 
The “ingoing” contacts on the upper hemisphere of the 
grain each carry a vertical force, which is analogous to an 
atom’s kinetic energy, and the collision within the grain 
redistributes these quanta into the “outgoing” contacts 
on the lower hemisphere. Thus, the set of all atoms in one 
layer of the q model represents the same set of atoms in 
the next layer, after they have emerged from their colli- 
sions. An TV x N lattice represents only N atoms, not N 2 . 
In summary, one spatial dimension becomes analogous to 
time while the other(s) represent the extent of the system 
at any one moment, so the set of just N atoms (contacts) 
within that instantaneous extent are propagated across 
the lattice in the remaining “time” direction. This anal- 
ogy can also be extended beyond the regular lattice of 
the q model to the arbitrarily disordered “lattice” of a 
real granular packing as illustrated in Fig. 2 if we relax 
the constraint that particle count is conserved during the 
collisions. 

This picture by Coppersmith, et al, is the beginning 
of a second analogy to thermodynamics that is comple- 
mentary to the one put forward by Edwards [5, 6]. In 
the thermodynamic analogy of Edwards, the geometry 
of the packing is made to explore its configuration space 
through a series of mechanical taps, and under certain 
conditions it eventually relaxes to the most entropic vol- 
ume, the one represented by the highest number of mi- 
crostates. Thus, the Edwards ergodicity is extrinsic and 
dynamic, affecting the contact geometry of the packing 
as well as its forces. The type of randomization that is 
pictured by the q model is intrinsic and static, affecting 
only the forces which propagate through the fabric, but 
achieving thermalization whether or not the contact ge- 
ometry has been allowed to relax. The problem with this 



FIG. 2: Illustration of the analogy by Coppersmith, et al. [4], 
between the static granular media and thermodynamics, in 
which the vertical dimension is analogous to the time dimen- 
sion. 

new analogy, however, is that it treats only one of the 
spatial dimensions as analogous to the time dimension, 
and thus it breaks the inherent symmetries of granular 
ensembles. Recovering that symmetry is what leads to 
the full, self-ergodic theory. 


C. Ergodicity in Orthogonal Dimensions 

The q model only deals with the vertical component of 
force (/ 2 ), conserving its sum layer-by-layer in the ver- 
sion of the q model without gravity, and illustrating a 
type of randomization (or ergodicity) in one conjugate 
spatial dimension. However, in a granular packing with- 
out gravity, the sum of the forces perpendicular to the x 
axis and to the y axis are also conserved layer-by-layer or 
column- by-column. To explain this statement and make 
it rigorous we must first define a “layer” of grain contacts. 
Referring to Fig. (3), we: (1) draw a cross-sectional line 
(like any one of the dashed lines in the figure), (2) se- 
lect the set of grains intersected by that one line (shaded 
grains), and (3) include only those contacts on that set 
of grains (heavy dots) that (a) connect to other grains 
external to the shaded set and (b) are located on one 
side of the shaded set. Next, we decompose the forces 
on that set of contacts into Cartesian components par- 
allel and normal to the dashed line. Summing all the 
components in the normal direction produces the “total 
Cartesian load” in that plane. In the absence of grav- 
ity, the total Cartesian load is conserved with respect to 
translation of the dashed line, as illustrated in the figure 
by the parallel layer of shaded grains and heavy dots. If 
that were not so, then the grains between the two shaded 
layers would be accelerating. 

Generalizing the key insight of the q model, a granu- 
lar packing should therefore have internal force relation- 
ships between the subspaces that represent not just the 
layer-by-layer subspaces in the z direction, but also the 
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FIG. 3: Illustration how the kinetic analogy of Coppersmith, 
granular packings. 

column-by-column subspaces in the x direction and again 
in the y direction (for 3D packings). Each such set of 
relationships imply their own randomization in the cor- 
responding direction subject to their own conservation 
law. 

There are no equivalent subspace relationships in a 
classical thermal system. In an ideal gas, for example, 
it is extremely unlikely that the sum of the particle mo- 
menta or energies in a cross-sectional layer of the con- 
tainer would be the same as in any other parallel cross- 
sectional layer at any one moment of time. Furthermore, 
dynamics cause the gas to evolve ergodically so that it 
visits every possible state equally without regard for any 
subspace correlations. This is obviously not the case for 
static granular packings, where the relationships between 
neighboring layer-wise subspaces must exist merely be- 
cause the packing is static. These internal relationships 
between the subspaces add a significantly new set of con- 
straints onto the form of the overall density of states. 

The thermodynamic analogy to the q model identifies 
the spatial dimension perpendicular to the layer as the 
time dimension. To extend the q model analogy to the 
full set of Cartesian components in a 2D granular pack- 
ing, we may consider both spatial dimensions to be time 
dimensions (a “plane” of time), and both independently- 
conserved quantities (total Cartesian loads) to be ener- 
gies. In 3D systems, we have three “time” dimensions 
and three independently conserved “energies” . These two 
(three) conservation laws are orthogonal to each other, so 
that we do not have a simple uni-directional ergodicity in 
the time dimension as in thermodynamics. Instead, we 
have two (three) randomizations (or ergodicities) that are 
orthogonal to one another. 


et al. [4], with respect to the q model, may be adapted to 2D 


D. Exploring Phase Space Without Moving 

The key insight of the q model is that contact forces 
can explore phase space in a granular packing that is 
completely static: the exploration of phase space is not 
performed by the entire system as a whole through the 
time dimension; rather it is performed layer-by-layer in 
the spatial dimensions because of the static equilibrium 
relationships that exist between the layers. 

To emphasize how unusual this concept actually is, we 
consider a granular packing’s phase space dynamics, or 
rather its lack of phase space dynamics. Each granu- 
lar packing is plotted in phase space as a single point, 
and it remains perfectly immobile wherever it may be. 
We can see that there is no exploration of phase space 
at all. However, only certain, very special locations in 
phase space are valid for static granular packings. These 
locations have a special property in that the subspaces 
of the coordinates (corresponding to the individual lay- 
ers of the packing) have static equilibrium relationships 
with one another. These relationships are such that each 
successive subspace has a locus that moves relative to the 
previous subspace according to the micromechanics of the 
layer of grains between them. These micromechanics con- 
serve the sum of the subspace coordinates in each layer 
while also introducing a randomization from one layer to 
the next. Thus, there is a random walk built right into 
every valid locus of phase space. It is not moving in time, 
but it really is exploring phase space, nonetheless. Be- 
fore a packing can become static it must first locate and 
settle upon one of these special loci that has the equiv- 
alence of motion built into itself, and not just in one 
dimension, but in two (or three) orthogonal decomposi- 
tions simultaneously. We shall further show below that 
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this orthogonal, ingenerate exploration of phase space is 
not just randomized, but really ergodic in the Boltzman- 
nian sense. This new type of ergodicity is a very impor- 
tant concept to understanding granular media and so it 
deserves a name. “Self-ergodicity” seems appropriate. 


E. Vectorized q Model 

Because these orthogonal randomizations (or self- 
ergodicity) appear to be the essential concept in solving 
granular contact force physics, it would seem important 
to create a version of the q model that treats the or- 
thogonal directions correctly. However, generalizing the 
q model to include horizontal forces presented an inter- 
esting challenge. In one example [7], the q model was 
vectorized so that the full vectorial forces are propagated 
layer-by-layer in one direction across the lattice from top 
to bottom, just as in the scalar q model. As the infor- 
mation is propagated, vectorial force balance is enforced 
on each grain. This conserves the sum of the vertical 
forces in each layer, but the horizontal forces grow with- 
out bounds as the information progresses through the lat- 
tice. This is partly because of the direction that the hor- 
izontal forces were being propagated across the lattice. 
The sum of the horizontal forces in a single column can- 
not be enforced unless the entire column is propagated 
column-by-column in the horizontal direction. Propagat- 
ing in that direction should allow the forces to explore 
phase space both up and down the column without bias, 
and this will statistically distribute them without bias 
in the upper versus lower parts of the lattice. This will 
solve the problem of mono tonic growth from top to bot- 
tom. Stated in its essence, granular ensembles in nature 
are fundamentally organized to conserve horizontal forces 
in the horizontal direction, and so to be fundamentally 
meaningful a model should also be organized this way. 

F. Overlaid q Models 

This observation indicates another way to attempt gen- 
eralizing the q model, illustrated in Fig. 4. For a 2D 
granular packing, we may simply implement two sepa- 
rate q models independently, throw away the upper, un- 
relaxed layers while keeping enough relaxed layers of each 
q model so that their lattices are square, and then rotate 
one of them by 7 t/ 2 radians with respect to the other. 
The two square lattices can then be sandwiched together 
so that each leg of the lattice is provided with two scalar 
forces, one the vertical component of force from the first 
q model, and the other the horizontal component of force 
from the rotated q model. This successfully achieves vec- 
torial force balance on each grain and it also organizes 
the global conservation laws properly in each orthogonal 
direction so that neither component of force grows with- 
out bounds. Thus, it solves the problem of the vectorial 
q model described above. 



FIG. 4: Two square, relaxed sections of q models may be over- 
laid to produce force vectors that (1) sum to zero on every 
grain, and (2) properly conserve the total force in each conju- 
gate direction. However, the correlations between force vector 
angles and contact surface angles are not properly organized 
in this model as illustrated in Fig. 5. 


Unfortunatly, that is not the end of the story. The 
resulting set of contact force vectors produced by this 
method is so wrong that it cannot even be considered 
a model of granular materials. In a real granular mate- 
rial, the directions of contact force vectors are correlated 
to the normal vectors at the contacting surfaces as il- 
lustrated in Fig. 5. For frictionless, cohesionless grains, 
the forces must lie perfectly on the contact normals. For 
Coulomb friction, the allowable force vector directions 
widen to become cones centered on the contact normals. 
For perfect friction but still no cohesion, the cones open 
further to become entire hemispheres looking inwardly 
toward the grain. For imperfect cohesion and perfect 
friction, the allowable directions become fully spherical 
for weak forces but still hemispheric for strong forces that 
exceed the allowable cohesion. Only in the case of perfect 
cohesion and perfect friction do the allowable directions 
of the contact force vectors become completely uncorre- 
lated to the contact normal directions, but that is pre- 
cisely the case where it is no longer a granular material. 
In general, the less fragile the packing, the less correlated 
the contact force directions are with the contact angles. 
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FIG. 5: In a realistic granular packing, the directions of con- 
tact force vectors cannot be random, but must be correlated 
to the contact surface normal directions. 


Fragility is what makes a granular material essentially 
granular [8]. In this overlaid q model the horizontal and 
vertical components of force are not correlated with each 
other to point in any special directions, and so the re- 
sulting network of forces cannot be interpreted in any 
meaningful way as a granular contact geometry. In fact, 
solving for grain diameters on the basis of this set of force 
vectors and assuming no cohesion, some of the grains are 
found to have negative diameters, clearly demonstrating 
that it is not geometrically meaningful. This method 
also produces the incorrect form for P/(/), beginning 
from P/( 0) = 0 rather than P/(0) > 0 as seen in the 
numerical simulations. This, too, is due to the unnatural 
statistical independence of the two force components at 
each contact. 


G. Competing Constraints in a Granular Packing 

Here is an indispensible insight. The key problem in 
generalizing kinetic theory to contact forces inside static 
granular materials is that the ergodicities are orthogonal 
but not independent. The orthogonal conservation laws 
require that the components of force information prop- 
agate in their own corresponding spatial directions, but 
both components of force must arrive at each grain simul- 
taneously in order to effect proper correlations between 
them at the contact angle. These two constraints are 
severe. Nature meets them, of course, by exploring the 
phase space through the real time dimension, in which 
the geometry of the grains’ packing adjusts itself until it 
finds one of the very rare stable states in force and geom- 
etry space that meet all orthogonal conservation laws in 
every layer and column and all the force vector direction 
constraints at every grain contact. 

Similarly, several lattice-based models have been de- 
vised to use the equivalence of temporal exploration of 
the phase space to find valid states. These models use 
either an iteration process on the lattice [9], or a global 


annealing process that leverages hyperstaticity [10], or an 
iteration that is localized to the first coordination shells 
by leveraging hyperstaticity [11]. Following specific itera- 
tions of a model in that way is analogous to following the 
state of a thermodynamic system through time, which is 
the idea in kinetic theory. The intention of this paper, 
on the other hand, is to derive a solvable Boltzmann-like 
equation and thus to obtain a statistical mechanics the- 
ory, not a kinetic theory. Therefore, it is essential to keep 
the specific iterations of the actual time dimension out 
of this theory. 

To develop a lattice model that avoids temporal iter- 
ations, we would need to propagate the waves of infor- 
mation in each orthogonal direction (for the conservation 
laws), and yet have each wave arrive at every grain in the 
packing simultaneously (to effect the force component 
correlations at the contacts). Of course, this is geometri- 
cally impossible; only a set of grains running diagonally 
across the lattice can receive both waves simultaneously. 
Therefore, we can conclude that it is impossible to gener- 
alize a q model in a way that is organized fundamentally 
the same as nature’s granular packings (without resort- 
ing to nature’s temporal iterations). The solution to this 
problem is to dispense with the lattice. We must col- 
lapse the density of multi particle states to a density of 
single particle states. That way, we can propagate the 
information to every particle in the density of states si- 
multaneously while enforcing all of the conservation laws. 
Instead of following Liouville or Gibbs with multi- particle 
phase spaces, we shall follow Boltzmann with single par- 
ticle distributions. This Boltzmannian approach allows 
us to “collide” grains together without knowledge of their 
spatial arrangement in the lattice. However, without any 
knowledge of the spatial arrangement, we will need the 
equivalent of Boltzmann’s stosszahlansatz to tell us how 
to deal with any spatial correlations that may exist. Since 
granular packings have the analog of multiple “time” di- 
mensions, this stosszahlansatz must be generalized in a 
way that properly maintains the extra symmetries of the 
physics. 


III. THE GENERALIZED STOSSZAHLANSATZ 

It will be shown below that the properly symmetric 
version of the stosszahlansatz is obtained automatically 
by taking the First Shell Approximation in the fabric in 
the density of multi particle states before collapsing it to 
a density of single particle states. We will take four steps 
to develop this idea. First, the density of multi particle 
states will be described from first principles as our start- 
ing point. Second, to understand how we must modify 
this density of states, qualitative arguments will be pro- 
vided to show that the First Shell Approximation does 
not discard essential force correlations in the packing, 
just as Boltzmann’s stosszahlansatz does not discard es- 
sential momentum correlations in the state of a dilute 
gas. In thermodynamics, Boltzmann’s stosszahlansatz 
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seems eminently reasonable, but a rigorous proof is lack- 
ing and so empirical testing is required. Likewise, the 
qualitative arguments provided here indicate that the 
First Shell Approximation is reasonable, but in the end 
testing will be required. Third, we will mathematically 
analyze the First Shell Approximation and finish with 
an expression that is recognizable by its form as the 
generalized version of Boltzmann’s stosszahlansatz. This 
will demonstrate that the First Shell Approximation re- 
ally does perform the same function for self-ergodicity 
as Boltmann’s assumption does for thermal ergodicity. 
Fourth, just as Boltzmann’s assumption in thermal sys- 
tems required empirical validation, so empirical valida- 
tion will be required for the generalized version. These 
shall be kept to the end of the paper after the theory’s 
predictions are fully developed. These results shall in- 
dicate that the generalized stosszahlansatz is excellent, 
perhaps to the same degree of rigor if not better than 
Boltmann’s original version in thermal systems. 

A. The Density of Multi Particle States 


to be isostatic in numerical simulations [16]. In this equa- 
tion, the geometrical information of the packing is con- 
tained in the set of variables £. The contact force on grain 
fi on its contact v is The construction list C tells 
which (/zi/) = (ary) connects to (pu/) = (/?<$), and hence 
is a part of the geometric information otherwise included 
in £. C is written explicitly because it will be impor- 
tant to the derivations later in this paper. The func- 
tion Q 40 computes the generalized fabric of the packing, 
#3 j 04 ) [17]. Although this P40 nomenclature 
implies that every grain has exactly Z = 4 contacts, P40 
may be understood more generally as a statistical sum- 
mary of all the contact angle information of the grains 
within their first coordination shells, regardless of their 
Z . The function 0^, introduced by Edwards [5], enforces 
non-overlapping of the grains such that each packing is 
physically realizable. Here, 0^ serves this purpose only 
in the second and higher coordination shells, since P 40 is 
sufficient to perform that same function within the first 
coordination shell. 

The duress tensor Dij is defined as the extensive coun- 
terpart to the stress tensor, 


The microcanonical density of multi particle states 
has been written [12-14] for static granular packings of 
rigid, round, frictionless grains (without the symmetry- 
breaking effect of gravity) in the form, 


p{ forces, geometry} = ]^[<S(Newton’s Third Law) 

contacts 

x 0(No Overlapping Grains) 
x (Newton’s Second Law) 

grains 

x [| 0 (No Tensile Forces) 

contacts 

x ^(Specified Stress State) 
x S (Specified Fabric State) . (2) 


with the expressions, 
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(3) 


It was adapted from Edwards’ microcanonical density of 
states in the geometric information [5] and from Edwards’ 
and Grinev’s stress probability functional in the force in- 
formation [15]. It assumes Edwards’ flat measure in the 
geometry and isostacy in the forces. Packings in the per- 
fectly frictionless, perfectly rigid idealization do appear 


= (Area) x cr^ (for 2D packings) 

= (Volume) x 0 {j (for 3D packings). (4) 

That is, it is a rank two tensor that has the units of 
force x distance and represents the total of all forces 
operating across their branch vector lengths within the 
volume of the packing. In the thermodynamic analogy, 
all the spatial dimensions are like time dimensions and so 
the multi-dimensional “duration” of the stress in a region 
of grains is the volume of that region. Hence, duration 
multiplied by stress is duress, which explains the name. 
It is related to the strain energy of an elastic system, 
although in a perfectly rigid granular packing there is no 
elastic strain energy. 

As explained in [12], it is possible in the idealized Z = 
4 case to change the coordinates of p from (f^ u X) to 
(w£,Wy,Q where the Cartesian loads are defined as 

z & 

< = ^E/H c ° s <n, < = ^E/h sin<n. (5) 

U~\ V 

In the Z = 4 case, specifying (w XJ w y ,9 1 ,. . . ,# 4) permits 
all four contact forces to be solved trigonometrically, f n = 
fn{g) for 7 ] = 1,2,3, and 4. For Z / 4 either more or 
fewer state variables should be specified. 


B. Qualitative Arguments for the First Shell 
Approximation 


Eq. 3 specifies the generalized fabric P40 as part of the 
definition of the microcanonical ensemble. This function 
imposes steric exclusion in the first coordination shell, 
because it evaluates to zero whenever 6 { and Oj ( i ^ j) 
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0£ enforces steric exclusion in P AQ (0 1# 0 2 , 9 3 , 0 4 ) enforces 
the second coordination shell steric exclusion in the first 

and higher coordination shell 

FIG. 6: Illustration of coordination shells. The central grain 
is shaded dark. The first coordination shell is shaded lightly. 
The second coordination shell is dashed and unshaded. Sim- 
ply omitting 0< from the density of states is the First Shell 
Approximation. 




FIG. 8: Correlations between neighboring contact forces may 
arise two different ways. (Left) They arise through intra-grain 
stability requirements, resulting in a strong, two-point corre- 
lation between the forces at B and A. (Right) They arise 
through a series of intra-grain stability requirements work- 
ing grain-by-grain around the loops, resulting in weak higher- 
order correlations. The example here shows a four-point cor- 
relation through the loop, B to C to D to A, which will be 
much weaker than the two-point correlation directly from B 
to A. 


No closure of force loops 




FIG. 7: (Left) When 0^ is kept as part of the density of states, 
forces form loops and so “colliding” forces on the central grain 
may have pre-correlation, because they have seen parts of each 
other before. (Right) When 0< is omitted from the density 
of states, which is the First Shell Approximation, then force 
loops do not close and so colliding forces are not pre-correlated 
by the packing. The only correlation arises from the stability 
requirements of the central grain, itself. 


are within 7r/3 of one another (as an example for the 
case of monodisperse spheres). The purpose of 0^ is to 
enforce steric exclusion in the second and higher coordi- 
nation shells, as illustrated in Fig. 6. Thus, eliminating 
from Eq. 3 is the First Shell Approximation, keeping 
steric exclusion in the first coordination shell via P 4 # but 
ignoring all geometric constraints in the shells further 
away. Fig. 7 shows that eliminating 0^ assures that the 
packings having closed force loops outside the first coor- 
dination shell will comprise a set of zero measure, not af- 
fecting the statistics of the ensemble. Therefore, we may 
say that (essentially) all forces arriving at a grain will 


have arrived directly from infinity without ever interact- 
ing with one another. Thus it ensures that no correlation 
between colliding forces can arise in the packing, apart 
from the stability requirements of the central grain itself. 
In the First Shell Approximation, all correlation is intra- 
grain and propagates outward through the packing; no 
correlation arises non-locally in the lattice propagating 
inward to shape the statistics of the central grain. This 
explains how the First Shell Approximation is equivalent 
to the generalized Boltmann assumption. 

Fig. 8 shows that the correlations due to those force 
loops are negligible and hence can be reasonably dis- 
carded. Correlations between neighboring contacts on 
the same grain arise either through the grain itself or 
through the loops in the packing. A typical loop is four or 
more grains, so going the long way around a loop induces 
a weak third-order correlation, but going the short way 
between the two contacts (staying intra-grain) induces a 
very strong first-order correlation. The First Shell Ap- 
proximation discards the very weak four point correlation 
but keeps the very strong two-point correlation. Thus, it 
is equivalent to truncating the BBGKY heirarchy, which 
in thermal statistical mechanics is equivalent in effect to 
Boltzmann’s stoss zahlansatz. 

It has been shown [2] that contact forces on the same 
grain are strongly correlated with one another. There is 
anti-correlation for contacts closer together than roughly 
7 t/ 2 radians of angular separation, and positive correla- 
tion when the angular separation is greather than roughly 
7 t/ 2. The correlation continues to increase as the con- 
tacts are increasingly distant from one another but still 
on the same grain. These strong intra-grain relationships 
make sense due to the requirements of static equilibrium 
of the individual grains. Contacts on the same quadrant 
compete for a share of the same load and hence are anti- 
correlated. Contacts opposite one another transmit load 
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FIG. 9: Contacts that are approximately 7r/2 radians away 
from one another on the same grain are only weakly correlated 
as illustrated by the closed loop of four grains that allows any 
combination of weak and strong force chains to pass through 
it. If the angles were precisely 7r/2, then the four force chains 
in this figure would be completely independent. 



FIG. 10: Contacts that are highly correlated are close to n 
radians apart on the same grain [2]. Hence, a closed loop 
composed of highly correlated pairs of contacts can turn only 
very slowly and must pass through a very large number of 
grains. 


through one another and hence are correlated. Simplis- 
tically we could expect n/2 to be the crossover point of 
no correlation as illustrated in Fig. (9). This is impor- 
tant, because small loops of grains are formed through 
adjacent pairs of contacts that are closer to n/2 than n 
radians of separation, and hence there is minimal cor- 
relation in each two-point leg of a force loop. The net 
four point correlation, which we have already argued to 
be weak, is therefore even weaker than what we would 
expect in a thermal system. On the other hand, force 
loops that pass through grains with closer to n radians 
of separation will require a vastly larger number of grains 
in the loop to slowly turn through the full 27r radians to 
close the loop as illustrated in Fig. 10, and these iV-point 
correlations (N >> 4) will be extraordinarily weak. 

The total correlation between a pair of contacts on a 


grain must be the sum of information from all the loops 
in the packing that contain the grain in question. Due 
to the disorder of the packing and the large number of 
loops that contain the same grain, some correlating and 
some anti-correlating, it is expected that the contribu- 
tions from increasingly larger loops of grains will be in- 
creasingly decoherent and largely cancel one another. 

From these arguments there is good reason to assume 
a priori that the intra-grain contribution to the correla- 
tions is the dominant one. Thus, it appears reasonable 
that the first shell approximation does not discard essen- 
tial correlations from a granular packing. 


C. Mathematical Derivation from the First Shell 
Approximation 

These physical arguments are the basic idea behind us- 
ing the first shell approximation as the stosszahlansatz, 
and now this idea shall be developed mathematically. 
The first step in collapsing the multi-particle density of 
states to a single particle density of states is to sum p 
over all grain exchanges represented by the N\ possible 
construction lists C/t, 


N! 

p{r , C I P*o,D ij )=H I C k ,P 4 e,D i:i } 

k = 1 

(6) 

where the tilde on p indicates that this density is in an 
assembly space rather than a phase space. Each locus in 
an assembly space specifies a set of single grain configura- 
tions, and the density in this space tells how many valid 
packing states may be assembled by connecting together 
the grains from that set. All terms in Eq. 3 factor out 
from the sum except 0^ and the product over the delta 
functions enforcing Newton’s Third Law, 


p{r, < i p ie ,D t] } = j^eac i c k )X[6\r+f 0S )\ 
U=i c k ) 

( N / Z° \ Z° ^ 


* n 

*,J=X.y \c*,0 y 

x <5(P 4 0-Q 4 »(O)- 


( 7 ) 


We call the summation in the brackets 3>, 


<*> = £ e< (c i c k ) n (n + P s ) (8) 

k= 1 C k 

We wish to determine the behavior of in the thermo- 
dynamic limit (very large N). Unless we discard Oq we 
can make no progress. All of the subsequent develop- 
ment results from this First Shell Approximation. We 
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may expand the product over C*, 

N\ N 4 N 4 

* - £ n n n n *(1* - t x ) <»> 

fc=l Q=1 7=1 0=0. <5=1 

where the Kronecker delta function Sc k is nonzero only 
for contact pairs contained within list C*. Note that the 
product in 0 is only from a to N, the upper triangle in 
the (a, 0) plane. This prevents enforcement of Newton’s 
Third Law twice for every contact, once looking from 
grain a to grain 0 and again from 0 to a, and produces 
the correct number of delta functions so that the units of 
p are correct. However, we actually want the symmetry of 
looking both ways since the ensemble is symmetric with 
respect to rotating the granular packings by n radians, 
and so in the next step we sum over the entire (a, 0) 
plane, doubling the number of delta functions, and we 
must write it as 4> 2 , 

N\ N 4 N 4 

* ! =,&Ennnn^(/ i " -/’")■ <><» 

fc=l a=l 7=1 0=1 6=1 


Next we commute the sum in k inside the two products. 


<F 2 


N 4 N\ N 4 


j-*<n nsnrK ‘H/- 

Q =1 7 =1 k=l 0=1 6=1 


Ot = 1 7=1 0=\ 6=1 



This introduces a need for renormalization with M be- 
cause each of the terms in the product is now a sum of 
many delta functions instead of just one delta function, 
vastly increasing the number of states in p. However, 
M will be suppressed henceforth since the normalization 
has no effect on the derivation. This commutation also 
introduces an error because the product of all these sums 
expands to a sum of the product of all possible combi- 
nations of the delta functions with replacement, whereas 
it should have been without replacement. In the limit 
N — > oo this error vanishes, because the majority of 
the terms in the product of delta functions will have the 
same representative distribution of grains p g , whereas the 
cross- terms with significantly deviating p g + Ap g will be 
a vanishingly small fraction and make zero contribution 
to the sum in that limit. Therefore, the renormalization 
of the increased number of states is trivial. 

Next, we expand the delta functions as the difference 
between two step functions with the limit of vanish- 
ing distance between them and with magnitudes propor- 
tional to the inverse of that distance. 

N 4 N 4 

$ 2 = hm nnn Hm i ^ 

N—>oc 11 It A/— *0 A / A0^O AO 

Q = 1 7 = 1 0 =\ 6=1 J 

* {e [f 0S - fi(D] - e [f» - fi+AD] } 

x {e[|0**-0 i (0°'»)|-ir] 

-©[I** 4 — 0 j+1 (<^)|— *]} (ii) 


0 




Integrate across 
time and point all 
momentum vectors 
inward 


7,+/,+/,+/.= o 


Identical to forces on 
static grains 


FIG. 11: Illustration how momentum vectors in binary col- 
lisions of a dilute gas are analogous to contact forces on a 
static grain. Boltzmann wrote the stosszahlansatz with pref- 
erence for the in-going momentum vectors, but for static gran- 
ular materials the symmetries among the contact force vectors 
must be restored. 


where the subscripts i and j define the i th interval of 
the force axis as the interval containing f ai and the j th 
interval of the angle axis as the interval containing 0 ai . 
Taking the summations, 


1 i n 4 

$ 2 = lim lim lim — — TT TT n X7 

N-+oo A f->0 A0— >0 AF A 6 3 


Q= 1 7=1 


( 12 ) 


where n X j is the number of grains in the summation that 
fall into the space between the step functions in the i th 
bin of the force axis and the j th bin of the angle axis. 
Taking the limits converts and taking the 

square root of both sides recovers $ in a usable form, 

N 

* = n T « ( i3 ) 

a=l 


where 


Ta = n /> / e /2 (/ a "^“ 7 ) (14) 

7=1 

or, simplifying notation, 


T “ = ri p }e\r) (15) 

7=1 

Now we are in a position to see that this is the properly 
symmetric generalization of Boltzmann’s stosszahlansatz. 
Referring to Fig. 11, in a particle collision the sum of the 
momentum vectors is conserved. Hence, if the outgoing 
momentum vectors are reversed, then the sum of all four 
momentum vectors will equal precisely zero. Putting a 
little space between the four vectors and drawing a circle, 
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we see that this set of four vectors is identical to four con- 
tact forces on a static grain. However, Boltzmann wrote 
the stosszahlansatz with preference for the forward time 
direction, and therefore he wrote that only the incoming 
pair of vectors are statistically uncorrelated. 

F(Pui>2,t) = f(pi,t) - f(p2,t) (16) 


If we try to apply this to granular contact forces,we could 
group the four forces on the grain into six different pairs, 
and there is no reason why one of those pairs should be 
written as uncorrelated to the exclusion of the other five. 
Comparing Eq. 15 to Eq. 16, we see that the former is 
a product over all four vectors, preserving the symme- 
try. Furthermore, the 1/2 power represents the analog 
of time-reversal symmetry. This is because the product 
in a includes all grains, so that a grain and its contact- 
ing neighbor are both included somewhere in the prod- 
uct. Thus, 1/2 power of a term is included with grain a, 
whereas the other 1/2 power of the same term is included 
with the contacting grain where the same force vector is 
looking in the opposite direction (time reversal symme- 
try). Thus, the generalized version of the stosszahlansatz 
not only treats all force vectors on a grain symmetrically, 
it also maintains the symmetry of Newton’s third law be- 
tween grains. 

For comparison with Eq. 15, consider another form 
that does not have the 1/2 power, 

T*(/i, f 2 , k, k) = P(fi) ■ Pih) ■ P(h) ■ P(k) (17) 


Why is T* not the correct generalization of Boltzmann’s 
stosszahlansatz ? Although it does maintain symmetry 
between all the contacts, it cannot be correct because it 
implies that all four forces are statistically independent. 
They cannot be, because static equilibrium has removed 
two degrees of freedom. On the other hand, Eq. 15 does 
not imply statistical independence because the 1/2 power 
ruins normalization and therefore P 1 ^ 2 (/) cannot be in- 
terpreted as a probability distribution. Furthermore, the 
removal of the two degrees of freedom through New- 
ton’s second law implies that the probability T of four 
specific forces colliding in a grain should have units of 
(Newtons) -2 rather than (Newtons) -4 as we see in T*. 
The 1/2 power in Eq. 15 solves this problem, too. 

In summary, we began this section with arguments 
that force correlations should not arise very strongly, if 
at all, through the closure of grain loops. This produced 
an expression T that may be interpreted as the probabil- 
ity that a set of four specific contact forces may “collide” 
with one another in a stable grain. This result will en- 
able the ergodicity proof for contact forces in the next 
section. 


TABLE I: Phase space conventions. 


Scale 

Number of 
Grains 

State 

Variable 


Density 
of States 

Grain 

1 

9 = (w x ,w y ,9i, 

• •■A) 

P(9) 

Set of Grains 

m 

II 

S 

5 9m) 

pin) 

Packing 

N 

F = ($1,52,- 

9N | C) 

p(0 


IV. THE SELF-ERGODIC THEOREM 
A. The Granular Transport Equation 

To simplify notation, Table I defines the conventions 
that describe the states of (1) single grains, (2) sets of 
grains, and (3) entire packings. The single-grain state 
variable g shown here assumes Z = 4 for every grain. 
This may be easily generalized for grains with Z ^ 4 but 
doing so adds no insight to the physics at this stage. 

Using these conventions, we proceed to develop a gran- 
ular contact force transport equation. Boltzmann wrote 
a transport equation to track the changes in the single 
particle distribution function f(v) through time. In writ- 
ing the analogous equation for contact forces we must use 
changes through space, not time, because self-ergodicity 
operates across space. We will therefore perform a trans- 
port of forces similar to the transport contained in the 
q model, following successive cross-sectional layers of 
grains. We use the word “layer” as it was defined with 
Fig. 3, so that “translating” refers to the continuous mo- 
tion of the cross-sectional line that picks out the set of 
grains in a layer. If the container is 2D, rectangular, 
and frictionless, then we will have two spatial axes along 
which to propagate layers of grains subject to the con- 
servation laws. We cannot perform the transport in a 
diagonal direction because then the length of layers will 
not be constant and neither will the total perpendicular 
force contained within them. This transport in a finite 
container is the microcanonical version of the transport 
because the system is closed. Alternatively, we can per- 
form the transport in an extremely large packing (per- 
haps infinite) and track only a finite cross sectional seg- 
ment that propagates across the packing in any random 
spatial direction, not only in the x or y direction. In 
that case, the conservation law will not be exact because 
some force will be constantly entering and exiting the 
segment at both of its ends. However, the longer the 
segment of grains is, then the less significant those fluc- 
tuations become and the more closely it approximates an 
exact conservation law. This is the canonical version of 
the transport, because the segment is in contact with an 
infinite reservoir of forces at each end. It does not mat- 
ter which way we imagine the transport to occur, nor in 
what direction. This will be re-addressed at the end. 

Suppose Maxwell’s demon chooses the layer of grains 
that we will use as our starting point. With malicious in- 
tent, he searches through the entire ensemble of packings 
and chooses a packing that has a single layer with a den- 
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sity of single grain states po(g) dissimilar to the density 
of its overall packing. That is, po(g) # p(</)- From this 
starting point, we wish to discover what happens to the 
density of layer states as a function of x as we propagate 
the cross sectional line in the x direction ( x being perpen- 
dicular to the layer), and so we write p^g) = p^g, x). As 
we translate from x — * x 4- Ax with Ax <<<< D par ticie 
(the grain diameter), a small fraction of the grains in 
the layer will no longer be intersected by the line and 
hence will exit the layer. Also, some new grains will be 
intersected by the cross-sectional line and hence they will 
join the layer. If the layer contains M grains, then the 
number of grains expected to leave the layer in Ax is 
m — MAx/Dparticle* If the fabric is constant across the 
packing (so that M is the same for every layer to good 
approximation) , then this is also the number of grains ex- 
pected to join the layer in Ax. For sufficiently small Ax 
the grains exiting the layer will be sufficiently far apart 
to be statistically independent. The probability for a 
particular set of m grains to leave during Ax is therefore 

m 

^out(T) = n Pl (5a) (18) 

a=l 

The probability for a particular set of m grains to en- 
ter during Ax can be written in terms of the general- 
ized stosszahlansatz from Eq. 15, except that we must be 
careful because in general P(f) for contacts that are on 
one hemisphere of the grains will not be the same as for 
the other hemisphere. That is because p(g) is evolving 
layer-by-layer, and P = P[p\- Therefore, we write the 
stosszahlansatz in the form, 

rHg) = n VP'(fv(9)) n VP(fv(g)) (w) 

77+ 7 ]- 

where rf~ refers to all the contacts on grain a in the 
reverse direction of the transport, and 77 + refers to 
its contacts in the forward direction of the transport. 
Likewise, P represents the distribution of external forces 
in the reverse direction and P f the distribution of con- 
tacts in the forward direction. P' is not known because 
it looks forward one step beyond the present state in the 
transport process. This difficulty is due to the inherent 
nature of granular materials: every grain touches a net- 
work containing every other grain, so there is no way 
to treat the transport process as causal in one direction, 
only. Fortunately, we do not need to know the form of 
P l to obtain the ergoic proof. Because the grains enter- 
ing the layer are statistically independent for sufficiently 
small Ax, we may write the probability for a particular 
set of m grains to enter during that interval as 

m 

^in(7')=^II T± (^) ( 20 ) 

Q = 1 

where N is for normalization. 

Because every new layer must exactly obey the con- 
servation law no matter how large or small Ax, we may 


consider this “in-out” structure to be analogous to Boltz- 
mann’s binary collisions, except that they are “m-nary” 
collisions and that each iteration Ax contains only one 
of them. The probability of having a particular collision 
7 — » 7 ' is 

771 

n ( 7-»v) = 

Q = 1 rf+ 

*T[>/P(fM) 

T)~ 

m 

= V(7) JI piga^ig'a) (21) 

a = 1 

where J\f is for normalization. The number of times that 
7 will go to any set of single grain states is 

- rri 

n ( 7 - all) = Vg[..- Vg' m V( 7 ) II P(^) T± (Sa) 

•'O Q=1 

(22) 

where the integrals are carried out only over the “stable” 
region © in &, that is, the region where none of the 
contact forces are tensile. Likewise, 

- TTL 

n(all — * 7 ) = / Vg\ -Vg' m ^( 7 ') H p(Sa) T± (Sa) 

a=l 

(23) 

Because every layer must connect to another layer, we 
insist that 

m 

n { 7 -» all) = H p(g Q ) (24) 

tt = 1 

and therefore 

r ~ 7n 

M~ l ( 7) = / Vg[ - Vg' m II^Q?;) (25) 

so that 

V(7)=A r[P,P')=Mtf) (26) 

which shall be important below. 

To determine the rate of change in p(g) during the 
transport process, we write, 

^P(5i) = J Vg 2 ■ ■ ■ Vg m [n(all -> 7) - M 7 -> all)) 

(27) 

B. Counting States for Granular Packings 

To evaluate how this transport equation behaves, we 
need a metric, similar to Boltzmann’s H , that will in- 
dicate how many packing states correspond to any par- 
ticular p(g). The “most entropic” p(g) is the one that 
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arises in the greatest number of packings. As explained 
in Ref. [12], the generalized stosszahlansatz enables the 
explicit counting of states as a functional of p(g). First 
we discretize p(g) — > i 'ijkimn, where the six arguments 
of g have been broken into small bins of size A w and 
A0 and indexed as flit, . . . ,# 4 n* Each bin is 

further divided into S smaller bins to enable the typi- 
cal binomial counting method for particles that do not 
have Pauli exclusion. S shall drop out of the equa- 
tions in the continuum and thermodynamic limits when 
A g — (Aw) 2 (A9) a — * 0 as N — ► oo. For now, with S 
we may write the number of states corresponding to a 
particuler v iikimn , 


nnnnrin 

i j k l m n 


' (S-1 + l/i...n)! ' 
.(5-1)! 


X 



0/ 

.71 * I...T1 


(28) 


Recall T was derived in Eq. 15 from Newton’s third law 
in Eq. 3 and is the fraction of grains at any locus in 
phase space that satisfy Newton’s third law with their 
neighbors. Likewise ’L was derived from the Newton’s 
second law in Eq. 3. It evaluates to zero at a locus of 
phase space where any of the g a cannot be stable apart 
from tensile forces, and it evaluates to unity elsewhere. 
We will divide the range of g, call it 5-space, into the 
two regions 


/ Vg= [ Vg+ f Vg (29) 

Jg Jo Jo 

where O contains all the grain configurations that have 
compressive contact forces and 0 contains all the grain 
configurations that have one or more tensile contact 
forces. Since we are dealing with cohesionless grains, O 
is the stable region and 0 is the unstable region. For 
simplicity we will restrict all further mathematics to 0 
so we may drop from the expressions. 

We define H as 


H± lim log ft (30) 

N —> oo 
A<?->0 

The natural logarithm of ft, using Sterling’s approxima- 
tion, may be written 

log Q = [(5 - 1 + Vi...n) log(S - 1 + Vi ... „) 

-(5-1) log (5 - 1) - Vi... n log Vi... n 
log * tj.] (31) 

Expanding the first term in a Taylor series around */*... n = 
0, setting S = NAg , and taking the continuum and ther- 
modynamic limits such that S » Vi... n , we obtain 

H = -fvgp(g) log ( 32 ) 


The functional H = H[p(g)\ is a measure of the number 
of states in p(F) that correspond to p(g ), and is in fact 
the generalization of Shannon’s entropy for granular con- 
tact forces. We may note that the H functional may be 
written, 

H — K - h (33) 

where 

h= [ Vg p(g)\ogp(g) (34) 

Jq 

looks more like Boltzmann’s original version of H , but 
where the correction term K accounts for the very con- 
straining, intimate connectedness of the grains in a pack- 
ing, 

K= [ Vg p(g) logT(g) (35) 

Jo 

Since T < 1, K < 0 and we see that Newton’s third law 
reduces the entropy of a granular packing as we would 
expect. 

During a transport process, p(g) must be allowed to 
evolve layer-by-layer, and therefore so must P(f). Hence, 
we distinguish between the hemispheres of a grain and 
use the form of T from Eq. 19 to write, 

H = -J°Vgp(g) log^y (36) 


C. Behavior of p and H in the Transport Process 


With this metric in hand, the question we wish to ad- 
dress is whether (d/d x)H > 0 during the transport pro- 
cess described above. Differentiating H , 



1 + log 


jHgY 
p(g) . 


+p(5)^ log r(s)} 


(37) 


The last term (call it x) may be expanded, 


X = J^Vg p{g)-^\ogT(g) 

= J^Vg p(g)^\ogl[i/P(Ug)) 

= \ J o v 9 Pig^ogPUM) (38) 


and then evaluated by changing the integration from a 
sum over all grains to a sum over all contacts, 


rOO J 

= i 

-S l d/P(/) 


(39) 
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Substituting Eqs. 22, 23, and 27 into (d/d x)H and then 
factoring out M as allowed by Eq. 26, 




- .Vg' m 


x 

x 



(40) 


Since we integrate over O for all g a and g ' a , we may swap 
variables and average the various equivalent expressions 
to obtain, 






X log 


riaPQ&Wga) 

n o P(gl) T± (ga) 


(41) 


By inspection of the integrand we see that it is never 
negative for any part of the region of integration. Hence, 


4-H > 0 (42) 

dx 


so that the entropy of the packing can only increase layer- 
by-layer. Furthermore, (d/d x)H = 0 if and only if 


n^a)T ± (fla) = n^) T± (fla) V 9a, 9a (43) 

a a 

and this is true if and only (d/d x)p(g) = 0 for all g. 
When that is the case, then P = P'. 

We can draw the following conclusions, which depend 
upon large M and the assumption that the generalized 
stosszahlansatz is valid. First, if there exists a layer in a 
granular packing that is not at maximum entropy, then 
the surrounding layers in both directions will (essentially 
always) relax to maximum entropy. However, beginning 
from any layer that is already relaxed at maximum en- 
tropy, then it is highly unlikely that its neighbors will 
fluctuate very far away from maximum entropy. Hence, 
the existence of any layer far from maximum entropy 
is exceedingly unlikely and such layers will practically 
never be found. Therefore, essentially all granular pack- 
ings must exist at maximum entropy. The propagation 
of contact forces in a granular material is truly ergodic 
in the special form that we are calling self-ergodic. 

This also proves that Edward’s flat measure in the ge- 
ometry of packings is not what determines the flat mea- 
sure in the contact forces. Suppose (for the sake of il- 
lustration) that Edwards’ flat measure were not correct 
such that some stable T states were more probable than 
others. The self-ergodicity theorem has proven that p(g) 
for those more probable F states must be identical to 


the p(g) for the less probable T states, as long as the 
stosszahlansatz is valid. Therefore, each packing will 
have the same ergodicity of contact forces regardless of 
its relative probability in T space. 

We may compare this to the Poincare Theorem, which 
tells us that there is a Gibbsian flat measure over the 
entire Poincare cycle (to within arbitrarily small dis- 
tances from every state). Whether there is or not such 
a Gibbsian flat measure in thermodynamics is unimpor- 
tant to Boltzmann’s proof. That proof tells us that any 
non-Maxwellian portion of the Poincare cycle is non- 
Maxwellian precisely because the stosszahlansatz was 
pathologically violated in the segment leading up to 
it. Because we know the stosszahlansatz is violated 
only rarely, not greatly, and not for long, this proves 
that essentially every part of the Poincare cycle has the 
Maxwell-Boltzmann distribution of velocities, whether or 
not the cycle itself has a flat measure. Thus, the dom- 
inance of the Maxwell-Boltzmann distribution is proven 
independently of Gibbsian flatness. Boltzmann’s practi- 
cal ergodicity does not depend upon ultimate mathemat- 
ical ergodicity. For exactly the same reason, the proof of 
self-ergodicity for granular contact forces and the deriva- 
tion of p(g) are independent of Edwards’ flat measure. 

However, this analogy is not exactly correct. Boltz- 
mann’s practical ergodicity tells us that every tiny 
segment of the Poincare cycle is dominated by the 
Maxwellian distribution, so that only a tiny amount of 
dynamics are needed to leave any non-Maxwellian state 
and settle into maximum entropy. So both the Poincare 
ergodicity and the Boltzmann ergodicity are attained by 
traveling along the trajectory in phase space. Static 
granular packings, on the other hand, do not travel 
along any trajectory in phase space. The self-ergodicity 
proven above depends upon contact forces maintaining 
static spatial relationships, not travelling and interacting 
through spatio-temporal ones. Thus, the self-ergodicity 
theorem does not depend upon any travelling through 
T space at all. Boltzmann’s proof obviated the need to 
travel the entire Poincare cycle and showed that the sys- 
tem need travel only a tiny segement of it. The self- 
ergodicity theorem on the other hand says that a granular 
packing need not travel through T space at all, because 
it has special relationships between its own subspaces 
built into its single locus in T space. This differentiates 
self-ergodicity from the ordinary ergodicity of thermody- 
namics. 


D. Derivation of the Density of States and P(f) 


The self-ergodicity theorem tells us that (d/d x)H = 0 
is the state of essentially every possible layer in every 
static granular packing and that the sufficient and neces- 
sary condition for this state is given by Eq. 43. This can 
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be written in the form, 


fr p(gq) n P&) 

iJi T± (ffo) 


v Sa.Sa 


(44) 


which implies that either side of the equation may be 
written as equal to a constant. Taking the logarithm, 


H>°g 

Q = 1 


Pij a) 

r*(9a) 


= c 


(45) 


and we see that is in the form of a conservation law. Its 
most general solution is when C is written as a linear 
combination of all conserved quantities. Since we did 
not specify the orientation of the layer or the direction 
of transport, this result is valid for all orientations and 
directions. Therefore, all components of the duress tensor 
from Eq. 4 must be conserved. 

Is fabric conserved in the transport process? By ob- 
vious geometry, each grain has a contact angle in the 
opposite direction of the corresponding contact in the 
adjacent layer. Thus, the distribution of contact angles 
in one layer P${6) must be related to the distribution in 
the adjacent layer Pq(0) by 


RAPe) = P'e (46) 

where R n is the operator for 7T-radian rotations. How- 

ever, it is not necessarily true that R^^P^e) = P±q. Thus, 
we cannot say a priori that 

K(p( 9)) = p'(g) (47) 

However, we can make another assumption additional to 

the stosszahlansatz. We assume that the construction 
method of the packing was such that the ensemble has 
Rn symmetry. Construction methods such as compact- 
ing, shearing, and shaking should satisfy that require- 
ment. Pouring in gravity does not. With this added 
assumption, (d/d x)p(g) = 0 implies that 


Rw(p4e) = P 49 (48) 

and 

Pie = Pie (49) 

so that fabric is conserved in the transport process. 
Therefore, the conserved quantities in the transport pro- 
cess are the same as the last two delta functions at the 
end of Eq. 3. 

Writing Eq. 45 in terms of the conserved quantities, 
we obtain for every g Q € O 

i° g + /W + ... A) (50) 

and solving for p, 

Pig) = G(0 J , . . . , e 4 )<b(g)r(g)e-^ (51) 


(Recall that # is the function that enforces the bounds on 
O, evaluating either to unity or zero if the cohesionless 
grain is stable or unstable, respectively.) G( 6 \, . . . , 64 ) 
is the fabric partition factor. Eq. 51 is identical to the 
equation derived in [12] by counting states and assuming 
maximum logfl, using Lagrange multipliers to conserve 
fabric and stress. Remembering that T is a functional of 
P(f) and that 

P(f) = f Vg p(g) 5 2 (f — f(g)) (52) 

Jo 

we see that Eq. 51 and 52 form a recursion so that p may 
be solved numerically when stress and fabric are spec- 
ified. Solving for p(g) provides everything that can be 
known about the density of single grain states, including 

PUY 


V. EMPIRICAL VALIDATION OF 
SELF-ERGODICITY 


The qualitative arguments of Sec. Ill told us that 
we could discard any force correlation possibly arising 
through force loops. This generalized stosszahlansatz 
must now be tested to show that any discarded corre- 
lations were truly inessential. We perform this test by 
solving for p(g) and comparing the results to numerical 
simulations of granular packings. 

For the present Eq. 51 has been solved in the isotropic 
case, assuming Z = 4 for every grain with the approxi- 
mation 


Pg (wx j Wy ,9p) p w (w x , vjy)pe (6p ) (w x , Wy )Q[ 3 ) (53) 

as described in [12]. This modified separability assumes 
no correlation between the loads and fabric apart from 
the truncating effect of 0$. The physical idea is that 
correlation does arise predominantly because nature dis- 
allows unstable grains. Empirical results have shown this 
to be correct [12]. In the remainder of this paper, “the 
theory” refers to the resulting numerical solution. 

The predictions of the theory have been overwhelm- 
ingly validated. A subset of these results have been re- 
ported earlier [18]. Fig. 12 shows P/{f) obtained from 
the theory and from the DEM data including all grains 
having Z > 2. The agreement is even better when only 
the Z = 4 grains from the DEM are included, as shown 
in Fig. 13. It should be noted that there are no free pa- 
rameters in the theory that could be adjusted to obtain 
these good fits. The predictions either fit the empirics or 
they do not. As it turns out, they are in such remarkable 
agreement that we may claim that the theory’s prediction 
is proven correct. 

Fig. 14 shows P x (fx) obtained from the theory and 
from the DEM data including all grains having Z > 2. 
Fig. 15 shows the comparison for the Z = 4 population of 
the DEM. Again, the prediction has been proven correct. 
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FIG. 12: Semi logarithmic Pj (/) obtained in numerical solu- 
tion of the transport equation and from DEM data, all Z > 2. 



FIG. 13: Semi logarithmic Pf (/) obtained in numerical solu- 
tion of the transport equation and from Z = 4 population of 
DEM. 


FIG. 14: Semi logarithmic P x ( f x ) obtained in numerical solu- 
tion of the transport equation and from DEM data, all Z > 2. 



FIG. 15: Semi logarithmic P x (f x ) obtained in numerical so- 
lution of the transport equation and from Z = 4 population 
of DEM data. 


It should also be noted that the forms for Poj(f,0) 
and P x (f x ) are conjugate to one another in that it is 
possible to convert back and forth between them using 
Youngquist’s transform (when P x (f x ) is known for all 
arbitrary rotations of the Cartesian axes, which for the 
isotropic case is trivial) [19]. Functional forms have been 
fitted to the empirical data for P/(f) and P x (f x ) from 
these simulations, and applying Youngquist’s transform 
does produce the correct functional forms to fit their con- 
jugate distributions. This transform proves analytically 
that when P/(f) has an exponentially decaying form, 
then P x (f x ) is a series sum of modified Bessel functions 
of the second kind. The two knees indicative of that type 
of Bessel function are clearly visible in Figs. 14 and 15. 

The shear ratio s = (w x — w y )/(w x + w y ) and the 
grain pressure t = w x + w y , calculated for each grain, are 
predicted by the theory to be statistically independent 
(unlike w x and w y ). The data from the DEM have been 
sorted into separate populations by Z so that, 

5 

Pst(s,t)=J2 n zP { s t\s,t) 

Z = 3 


where nz is the fraction of grains within each Z pop- 
ulation and Pst\ s it) is the distribution of each popu- 
lation taken separately. Then, we find that the s and 
t parameters are indeed statistically independent within 
each population considered separately, 

PsP(s,t) = P^\s)Pt Z \t) (55) 

for each value of Z. Again, the predictions of the theory 
have been validated. 

Fig. 16 compares the theory’s prediction and the DEM 
data for the distribution of the shear ratio P s {s). The 
agreement is quite good, even without sorting by Z. 
Since the theory assumed Z = 4 we wish to test whether 
the other Z populations fit the same form. For compar- 
ison we use the functional form that was fitted to the 
theory’s solution, 

P s (s) = cos(7rs/2) e _8s2 . (56) 

Remarkably, we find that the same form fits all three 
populations reasonably well after we write the standard 
deviation as 


(54) 


a = l/Z 


(57) 
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FIG. 16: P$(s) obtained in numerical solution of the transport 
equation and from DEM data including all grains having Z > 
2 . 



FIG. 17: P 8 (s) in DEM data segregated by Z. 


This form compared to all three Z populations is shown 
in Fig. 17. 

Examining the distribution of grain pressures, P t (i) 
where t = w x + w y , Fig. 18 compares the Z — 4 popu- 
lation of the DEM with the theory’s prediction. Fig. 19 
makes the same comparison semi logarithmically. These 
demonstrate good agreement. However, the DEM data 
with all Z populations included does not match the the- 
ory’s predictions very well due to the presence of the 
Z = 3 and Z = 5 grains. Nevertheless, as was done for 
P s (s) i an interesting comparison can be made between 
all these populations by finding a functional fit and pa- 
rameterizing it in terms of Z. A function that fits the 
theory’s prediction reasonably well over the region of ex- 
perimental interest (0 < t < 3), is 

P t (t) = tP-'e- 01 (58) 

with (3 = 5. After segregating the DEM data by Z we 
find that all three populations do fit this same form as 
shown in Figs. 20 and 21, but using (3 = 2Z — 4. 

This provides empirical proof that the generalized 
stosszahlansatz does not throw away any of the corre- 



FIG. 18: Pt(t) obtained from numerical solution of the trans- 
port equation and from DEM data including only grains hav- 
ing Z = 4. 



FIG. 19: Semi logarithmic comparison of P t (t) obtained from 
numerical solution of the transport equation and from DEM 
data including only grains having Z = 4. 

lations responsible for the major features of p g , just like 
Boltzmann’s original version does not throw away the es- 
sential features of the density of single particle states for 
the molecules in a dilute, thermalized gas. 

VI. SUMMARY AND CONCLUSIONS 

This paper has proven that the propagation of gran- 
ular contact forces is a truly ergodic process, although 
the nature of that ergodicity is more complex and in- 
teresting than what we find in thermal systems and so 
the term “self-ergodicity” appears to be a more appro- 
priate description. On the way to making this proof, it 
generalized Boltzmann’s stosszahlansatz and Shannon’s 
entropy for granular contact forces. As a consequence 
of the proof, it found a second way to derive p(g ), ob- 
taining agreement with the earlier method developed by 
the author. This p(g) is in outstanding agreement with 
empirical observations. 

Because the percolation of forces is self-ergodic, essen- 
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FIG. 20: Pt(t) in the DEM, segregated by Z. 
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tially every granular packing of the class considered in 
this paper will exist in a state of maximum entropy, with- 
out first having to explore phase space at all. There are 
of course many significant classes of granular packings 
that will not be relaxed to maximum entropy, but these 


have not been discussed here. They would include gran- 
ular packings that were prepared to have abrupt changes 
in fabric somewhere within their bulk so that the den- 
sity of states cannot be in its most relaxed state within 
the narrow band of grains on both sides of the interface. 
Another case would be the thin layer of grains at the 
top surface of a packing in gravity, where the self-weight 
introduces non-negligible stress gradients. However, so- 
lution of the more symmetric case considered here opens 
the door to solving more complex problems. 

While this paper has addressed frictionless, isostatic 
packings exclusively, early tests of frictional packings 
demonstrate that the theory will apply surprisingly well 
there as well, with only minor modifications to include 
the effects of the frictional forces. It also appears easy 
to modify the theory to variable coordination number 
[13, 18]. In brief, the theory appears to be rigorous and 
widely applicable and should form the basis to solve a 
wide array of problems in granular statistical mechanics. 

One statistical mechanics application in particular was 
first suggested by the author several years ago [20]. The 
idea is that, since granular contact force propagation is 
truly ergodic, then it is possible to define a contact force 
“temperature” and even a coordination number “chemi- 
cal potential” to explain the partition of stress and fabric 
fluctuations throughout a granular packing, and possibly 
leading to a full theory of rheology. These observations 
were apparent when numerical solutions to the theory 
first proved robust and convergent, even before this for- 
mal proof of ergodicity was accomplished. That is be- 
cause the numerical convergence discussed in the earlier 
publications was in fact an empirical proof of ergodic- 
ity, just as valid as the formal proof found in this paper. 
Once ergodicity is known, and once its characteristics are 
understood, then the temperature and chemical potential 
concepts fall out rather straightforwardly. A future pub- 
lication will be forthcoming to fully explain these con- 
cepts along with a series of discrete element modeling 
simulations that have been performed to test them and 
to draw further conclusions. 
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